* ============================================
$TITLE: Cournot, Pigou and Ricardo walk in a bar
* Unilateral environmental policy and leakage with market power and firm heterogeneity.
* Claudio Baccianti (Tilburg University) and Oliver Schenker (Frankfurt School of Finance & Management)

* This version: 2021-11-04
* ============================================

** This is the final version used in the paper
** Includes an homogenous sector

SCALAR vtol  tolerance level for convergence /1e-6/;

** options
SCALAR   PRINT   'save in Excel' /1/,
         FLEX_E  'emissions change with AA' /0/;


* -------------------------------------------
SETS     r       Countries /H, F/
         rw(r)   country with endogenous wage /F/,
         j       Sector /AGR, IND, MAN/,
         i       table item /GHG, GVA/;

ALIAS (r,rr,s);

SET oth(r,rr);
oth(r,rr) = YES$(ORD(r) NE ORD(rr));
* -------------------------------------------

* -------------------------------------------
* parameters calibrated from external sources

* Demand elasticities
PARAMETERS sigma(r)         Elasticity of substitution between varieties /H  4, F  4/,
           sigma0(r)        Elasticity of substitution between HG and DG /H  0.8, F  0.8/,
           alpha(r)         CES elasticity parameter differentiated varieties ,
           eta(r);

*** compute other parameters
alpha(r) = (sigma(r) - 1)/sigma(r) ;
eta(r)   = (1-alpha(r))/alpha(r);

* Pareto distribution and cost parameter of firms
PARAMETER k(r)     Shape param   /#r 4.5/,
          z0(r)    Lower bound interval pareto distr /#r 1/,
          delta    Death rate /0.02/,
* Cost parameters of firms
          gamma(r) Cost share of emissions /#r 0.01/,
* Total Factor Producitivity
          AA(r)    Aggregate productivity - differentiated goods,
          Ah(r)    Productivity - homogeneous good /#r 1/,
          A0(r)    Benchmark aggr productivity /#r 1/;

AA(r) = A0(r);
Ah(r) = 1;

* Climate policy
PARAMETER scc        'Social cost of carbon for Home' /95/,
          cp0(r)     'Initial carbon price (EUR/ton)'
          flagcp(r)  'Indicator 1/0 for carbon costs' /H 1, F 0/,
          flgcbam(r) 'Indicator for cbam' /H 0, F 0/;

* ================================================================
*** Data

PARAMETER tau(r)        'Variable trade costs for export market' /#r 1.04/,
          shH(r)        'Target expenditure share - Homogeneous goods',
          Ls(r)         'Labour supply (in million units)',
          gdp0(r)       'GDP in bn EUR'    /#r 13953/,
          SHEX(r)       'Benchmark share of exporters' /#r 0.5/,
          imgdp(r)      'Benchmark import penetration (on GDP)' /#r 0.435/,
          MKP0(r)       'Benchmark markup nontraded products' /#r 1.6/,
          prf_gdp0(r)   'Profits on gdp' /#r 0.31/,
          m0            'Varieties - Scaling factor',
* parameters calibrated from external sources
          n0(r)          'Number of firms',
          lambda(r)      Fixed costs - operations ,
          lambda_x(r)    Fixed costs - exports ,
* parameters to be calibrated targeting data
          eta_d          Coefficient of damage function (for Gt units),
          FE(r)          Entry fixed costs,
          pi(r)          Expenditure share parameter ,
* combination of parameters
          coeff_av(r)      combination of coefficients linking average to cutoff ,
          coeff_av2(r)     combination of coefficients,
          exp1(r)           ,
* Emissions of differentiated good sector taken from ../data/04_EEA_emissions_share_diff_goods.dta
          E0(r)            Benchmark level of sectoral emissions (bn mtons CO2eq) /#r 1.62/,
          Rtil0(r)         Benchmark aggregate costs for firms
;


* Using Rauch's (1999) goods classification and his shares
* He computes that differentiated goods represent between 64.6 and 67.1% of total US manufactures
shH(r) = 1-0.6585;

* and because of wage = numeraire,
Ls(r) = (1 - gamma(r)*(1 - shH(r)))*gdp0(r);

eta_d = scc/(2*SUM(r,E0(r)));
n0(r) = MKP0(r)*(1-alpha(r))/(MKP0(r) - 1);
m0 = 451.433;

coeff_av(r) = (k(r)*(1-alpha(r))/((1+k(r))*(1-alpha(r)) - 1));
coeff_av2(r) = (-k(r)*(1-alpha(r))/((1+k(r))*(1-alpha(r)) - 1));
exp1(r) = (alpha(r) - k(r)*(1-alpha(r)))/(1 - alpha(r));

DISPLAY alpha, AA, n0, eta_d, exp1, coeff_av2, Ls, z0, E0, shH;

ABORT$(k("H")*(1-alpha("H")) - alpha("H") < 0) "Unsuitable distribution, H", k;
ABORT$(k("F")*(1-alpha("F")) - alpha("F") < 0) "Unsuitable distribution, F", k;
ABORT$(n0("H") < 0) "Negative n0", n0;



* =================================================
* MODEL


VARIABLES
OBJ_VAR,
CTAX        Optimal carbon tax ;

POSITIVE VARIABLES
Rtil(r)      Revenues net of markups - average firm (bn EUR),
Ptop(r)      Aggregate price index (NEW),
P(r)         Price index - differentiated goods,
w(r)         Wage index,
M(r)         Share of active varieties,
n(r)         Number of firms per variety
E(r)         Aggregate emissions (Gt),
CTAX         Optimal carbon tax ,
z_star(r)    Productivity cutoffs - survival,
zx_star(r)   Productivity cutoffs - export,
theta(r)     inverse markup in domestic market,
theta_x(r)   inverse markup in export market  ,
beta(r)      Expenditure share of HG sector,
LH_SH(r)     Share of labour to the HG sector,

* additional variables for calibration
PAR_FE(r)    parameter - entry costs (-> profits on GDP),
PAR_cp0(r)   benchmark unit cost of carbon (-> E)
PAR_l(r)     parameter - fixed costs  (-> GDP)
PAR_lx(r)    parameter - export fixed costs (-> share of exporters)
PAR_PI(r)    parameter - share parameter HG sector
;

* -------------------------------------------------------

$macro theta_a(r)        {(n(r) - 1 + alpha(r))/n(r)}
$macro coeff_mkp(r)      {(n(r) + SUM(rr$oth(r,rr),n(rr)) + alpha(r) - 1)/n(r)}
$macro z_star_xr1(r)     {SUM(rr$oth(r,rr),zx_star(rr))}
$macro cprice(r)         {(PAR_cp0(r) + flagcp(r)*CTAX)}
$macro ratio1(r)         {coeff_av2(r)*(z_star(r)*z_star_xr1(r))**k(r)/(z_star_xr1(r)**k(r) - z_star(r)**k(r))* (z_star_xr1(r)**exp1(r) - z_star(r)**exp1(r))}
$macro zeta_a(r)         {(ratio1(r)/z_star(r)**((1-alpha(r))/alpha(r)) - 1)}
$macro zeta(r)           {alpha(r)/((1+k(r))*(1-alpha(r)) - 1)*(z0(r)/z_star_xr1(r))**k(r)}
$macro zeta_x(r)         {alpha(r)/((1+k(r))*(1-alpha(r)) - 1)*(z0(r)/zx_star(r))**k(r)}
$macro omega(r)          {(z_star(r)/z_star_xr1(r))**k(r)}
$macro zbarm(r)          {(ratio1(r))**((1-alpha(r))/alpha(r))}
$macro zbarc(r)          {sum(rr$oth(r,rr), zx_star(rr))*coeff_av(r)**((1-alpha(r))/alpha(r))}
$macro chi(r)            {(z_star(r)/zx_star(r))**k(r)}
* It is named agginc but it is only the income from the DG sector
$macro aggrinc(r)        {(w(r)*Ls(r) + cprice(r)*E(r)) + sum(rr$oth(r,rr), flgcbam(rr)*(cprice(r)-cprice(rr))*CBAMEDE(rr))}
$macro aggrexpdg(r)      {(w(r)*Ls(r) + cprice(r)*E(r))*(1 - beta(r))}
$macro incdg(r)          {(w(r)*(1 - LH_SH(r))*Ls(r) + cprice(r)*E(r))}
$macro c_tilde(r)        {(1/AA(r))*(cprice(r)/gamma(r))**gamma(r)*(w(r)/(1-gamma(r)))**(1-gamma(r))}
$macro costdisd(r)       {n(r)*c_tilde(r)/(SUM(rr$oth(r,rr), ttc(rr)*n(rr)*c_tilde(rr)) + n(r)*c_tilde(r))}
$macro costdisx(r)       {ttc(r)*n(r)*c_tilde(r)/(SUM(rr$oth(r,rr), n(rr)*c_tilde(rr)) + ttc(r)*n(r)*c_tilde(r))}

* Per-unit CBAM
* Computed as difference in carbon price per unit of production.
$macro CBAM(r)           {(sum(s$oth(r,s),cprice(s))-cprice(r))*gamma(r)*c_tilde(r)/cprice(r)}
* Total trade costs
$macro TTC(r)            {(tau(r)+flgcbam(r)*cbam(r))}
* Imported emissions from r that are covered by CBAM
$macro CBAMEDE(r)        {451.433*M(r)*chi(r)*n(r)*gamma(r)*theta_x(r)*w.l(r)*par_lx(r)*(coeff_av(r)**((1-alpha(r))/alpha(r)))**(sigma(r)-1)/(cprice(r)*(1-theta_x(r)))}

** Macros for output of results
$macro omega_out(r)      {(z_star.L(r)/SUM(rr$oth(r,rr),zx_star.L(rr)))**k(r)}
$macro ratio1_out(r)     {coeff_av2(r)*(z_star.L(r)*SUM(rr$oth(r,rr),zx_star.L(rr)))**k(r)/(SUM(rr$oth(r,rr),zx_star.L(rr))**k(r) - z_star.L(r)**k(r))* (SUM(rr$oth(r,rr),zx_star.L(rr))**exp1(r) - z_star.L(r)**exp1(r))}
$macro cprice2(r)        {(PAR_cp0.L(r) + flagcp(r)*CTAX.L)}
$macro cbam2(r)          {{(sum(s$oth(r,s),cprice2(s))-cprice2(r))*gamma(r)*{(1/AA(r))*(cprice2(r)/gamma(r))**gamma(r)*(w.l(r)/(1-gamma(r)))**(1-gamma(r))}/cprice2(r)}/cprice2(r)}
$macro TTCl(r)           {(tau(r)+flgcbam(r)*cbam2(r))}


* ===========================================================================
EQUATIONS
* model
OBJ              For NLP,
EQN_FE(r)        Free entry condition,
SS(r)            Steady state number of varieties,
EQN_R(r)         Equation R tilde,
EQN_E(r)         Emissions,
EQN_cutoffs(r)   Cutoff,
EQN_P(r)         Price Index Differentiated Good,
EQN_Ptop(r)      Price Index Top-Level,
FACTOR_CLR(r)    Clearing Factor Market,
MKP_D(r)         Domestic markup,
MKP_X(r)         Export markup,
EQ_Beta(r)       Share of income spent on homogenous good
MKT_HG           Market Clearance Homogenous Good,
TRADE_BAL        Trade balance
CONSTR(r)        Auxiliary Constraint,
CONSTR2(r)       Auxiliary Constraint,
OPT_TAX          Optimal emission tax,

* equations for calibration
CAL2(r)          profits on GDP ,
CAL4(r)          share of exporters
;
* ===========================================================================
OBJ..            OBJ_VAR =E= 0;

* Free entry. Govern n. Currently Eq 21 in paper
EQN_FE(r)..      delta*PAR_FE(r)*(z_star(r)/z0(r))**k(r) =E= PAR_l(r)*((1-omega(r))*zeta_a(r) + omega(r)*zeta(r))
                         + PAR_lx(r)*zeta_x(r);

* Steady state flow of entries. Governs M. Currently Eq 20 in paper
SS(r)..          delta*M(r) =E= (1-M(r))*(z0(r)/z_star(r))**k(r);

* Revenue
EQN_R(r)..       Rtil(r)/w(r) =E= (1 - omega(r))*(theta_a(r)/(1-theta_a(r)))*PAR_l(r)*ratio1(r)/z_star(r)**(alpha(r)/(1-alpha(r))) +  omega(r)*(theta(r)/(1-theta(r)))*PAR_l(r)*coeff_av(r) + chi(r)*(theta_x(r)/(1-theta_x(r)))*PAR_lx(r)*coeff_av(r);

* Emissions
EQN_E(r)..       cprice(r)*E(r)/(m0*M(r)*n(r)) =E= gamma(r)*Rtil(r);

* Productivity cutoff link
EQN_cutoffs(r)..     zx_star(r) =E=  ((PAR_lx(r)/PAR_l(r))*(aggrexpdg(r)/SUM(rr$oth(r,rr),aggrexpdg(rr)))*((1-theta_a(r))/(1-theta_x(r)))**2)**((1 - alpha(r))/alpha(r)) * (ttc(r)*theta_a(r))/theta_x(r)*(P(r)/SUM(rr$oth(r,rr),P(rr))) * z_star(r);

* Top-level price index
EQN_Ptop(r)..    Ptop(r) =E= (PAR_PI(r)**sigma0(r) + ((1-PAR_PI(r))**sigma0(r))*P(r)**(1-sigma0(r)))**(1/(1-sigma0(r)));

* Price index differentiated good
EQN_P(r)..       P(r) =E=  (m0*M(r))**(-eta(r)) * c_tilde(r)* ((1-omega(r))*(zbarm(r)*theta_a(r))**(1/eta(r))
                           +  omega(r)*(zbarc(r)*theta(r))**(1/eta(r)))**(-eta(r));

* Factor market clearing
FACTOR_CLR(r)..  incdg(r)/(m0*M(r)*n(r)*w(r)) =E=
                         (1 - omega(r))*(PAR_l(r)/(1-theta_a(r)))*ratio1(r)/z_star(r)**(alpha(r)/(1-alpha(r))) +
                         omega(r)*(PAR_l(r)/(1-theta(r)))*coeff_av(r) +
                         chi(r)*(PAR_lx(r)/(1-theta_x(r)))*coeff_av(r);

MKP_D(r)..        theta(r) =E= costdisd(r)*coeff_mkp(r) ;

MKP_X(r)..        theta_x(r) =E= costdisx(r)*coeff_mkp(r);

* Fraction of income spent on homogeneous good
EQ_beta(r)..      beta(r) =E= (PAR_PI(r)**sigma0(r))*(Ptop(r))**(sigma0(r) - 1);

MKT_HG..         SUM(r, Ah(r)*LH_SH(r)*Ls(r)) =E= SUM(r, beta(r)*aggrinc(r)) ;

* trade deficit in DG sector must be offset by a trade surplus in HG
TRADE_BAL..    -M("H")*m0*n("H")*w("H")*chi("H")*(PAR_lx("H")/(1-theta_x("H")))*coeff_av("H")
                         + m0*M("H")*n("F")*w("F")*omega("H")*(PAR_lx("F")/(1-theta_x("F")))*coeff_av("H")
                 =E= beta("H")*aggrinc("H") - Ah("H")*LH_SH("H")*Ls("H");

CONSTR(r)..      z_star(r) =G= z0(r);

CONSTR2(r)..     zx_star(r) =G= z_star(r);

OPT_TAX..        PAR_cp0("H") + CTAX =E= 2*eta_d*SUM(r, E(r)) * Ptop("H") ;

CAL2(r)..        prf_gdp0(r)*incdg(r) =E= w(r)*PAR_FE(r)*m0*(1-M(r))*n(r);

CAL4(r)..        shex(r) - (z_star(r)/zx_star(r))**k(r) =E= 0;

* ===================================================================
* Define CALIBRATION MODEL
MODEL ASYM_CAL /OBJ, SS, EQN_FE, EQN_R, EQN_E, EQN_cutoffs, EQN_Ptop, EQ_beta,
                 EQN_P, FACTOR_CLR, MKP_D, MKP_X, CONSTR, MKT_HG, TRADE_BAL,
                 CAL2, CAL4/;
* ===================================================================

** FIX VARIABLES
w.FX("F")  = Ah("F")/Ah("H");
w.FX("H")  = Ah("H");
CTAX.FX    = 0;
E.FX(r)    = E0(r);
beta.FX(r) = shH(r);

* Lambda is set to constant
PAR_l.fx(r) = 1;

n.l(r) = n0(r);

z_star.L(r) = z0(r)*1.5;
zx_star.L(r) = z0(r)*2;

* initialize other variables
P.L(r)       = 0.05;
Ptop.L(r)    = 1;
M.l(r)       =  ((z0(r)/z_star.l(r))**k(r))/(delta+(z0(r)/z_star.l(r))**k(r));display M.l;
M.up(r)      = 1;
Rtil.L(r)    = {(n.l(r) - 1 + alpha(r))/n.l(r)}*(1-beta.L(r))*gdp0(r)/(m0*M.L(r)*n.l(r));
theta.L(r)   = 0.85;
theta_x.L(r) = 0.93;

LH_SH.UP(r) = 1;
LH_SH.LO(r) = 0;
LH_SH.L(r) = shH(r);

z_star.LO(r) = z0(r);
zx_star.LO(r) = z0(r);

PAR_cp0.L(r) = gamma(r)*(1-beta.L(r))*gdp0(r)/E0(r);
PAR_lx.l(r) = PAR_l.l(r);
PAR_FE.l(r) =  prf_gdp0(r)*(1-LH_SH.L(r))*gdp0(r)/(m0*(1-M.L(r))*n.l(r));
PAR_PI.UP(r) = 1;
PAR_PI.L(r) = beta.L(r)**(1/sigma0(r))*(Ptop.L(r)**((1-sigma0(r))/sigma0(r)));

*ASYM_CAL.tolinfeas = 1e-3;

* ========================================================
* Calibration model
SOLVE ASYM_CAL MINIMIZING OBJ_VAR USING NLP;
abort$(ASYM_CAL.objval > vtol) "Benchmark is inconsistent.";

* =======================================================
* store results of calibration
cp0(r)      = PAR_cp0.L(r);
fe(r)       = PAR_FE.L(r);
lambda_x(r) = PAR_lx.L(r);
lambda(r)   = PAR_l.L(r);
pi(r)       = PAR_PI.L(r);

*recalibrated climate damage valuation

eta_d = scc/(2*SUM(r,E0(r)) * Ptop.l("H"))

parameter theta_a_check(r); theta_a_check(r) = (n.l(r) - 1 + alpha(r))/n.l(r);

** CHECKS
ABORT$(theta.L("H")> theta_x.L("H")) "Export fixed costs too low", theta_x.l;
ABORT$(theta.L("H")< theta_a_check("H")) "Problem w operating fixed costs traded", theta.l;

DISPLAY fe, cp0, eta_d;


* ========================================================
* Model definitions
* BENCHMARK MODEL WITHOUT UNILATERAL TAX
MODEL ASYM_BCM /OBJ, EQN_FE, SS, EQN_R, EQN_E, EQN_cutoffs,
                 EQ_beta, EQN_Ptop, MKT_HG, TRADE_BAL,
                 EQN_P, FACTOR_CLR, MKP_D, MKP_X, CONSTR/;

* UNILATERAL OPTIMAL EMISSION POLICY MODEL
MODEL ASYM_OPT /OBJ, EQN_FE, SS, EQN_R, EQN_E, EQN_cutoffs, EQN_P,
                 EQ_beta, EQN_Ptop, MKT_HG, TRADE_BAL,
                 FACTOR_CLR, MKP_D, MKP_X, CONSTR,OPT_TAX/;
* =======================================================

SET ITR_AA /ITR_AA1*ITR_AA9/;

PARAMETERS SCAA(ITR_AA), xscal, BCMRES(*,*,*), OPTRES(*,*,*), CBAMRES(*,*,*), CTIL(r), CDI_hd, ef, DCMP(*,*,*), CHANGE_OPT(*,*,*,*), CHANGE_CBAM(*,*,*,*),DECOMP(*,*);

PAR_FE.fx(r) = fe(r);
PAR_lx.fx(r) = lambda_x(r);

*$ONTEXT
LOOP(ITR_AA,

**       Set up new benchmark level
         xscal = 4 - ORD(ITR_AA);

         IF(xscal > 0,
                 AA("F") = A0("F")*(1 + 0.02*xscal);
         ELSE
                 AA("F") = A0("F")/(1 - xscal*0.02);
         );

*        Adjust F emissions following AA
         IF(FLEX_E,
                 E0("F") = E0("H")*AA("F");
         );

         SCAA(ITR_AA) = xscal;
         DISPLAY xscal,AA, LS;

**       Simulate new benchmark level
         CTAX.FX = 0;
         E.FX(r) = E0(r);
         beta.FX(r) = shH(r);

*        must keep cp0 flexible
         PAR_cp0.UP(r) = inf;
         PAR_cp0.lo(r) = 0;
         PAR_cp0.L(r) = gamma(r)*gdp0(r)/E0(r);

*        keep pi flexible
         PAR_PI.UP(r) = 1;
         PAR_PI.lo(r) = 0;
         PAR_PI.L(r) = beta.L(r)**(1/sigma0(r))*(Ptop.L(r)**((1-sigma0(r))/sigma0(r)));

* ---------------------------------------------------------------------
         SOLVE ASYM_BCM MINIMIZING OBJ_VAR USING NLP;
         abort$(ASYM_BCM.objval > vtol) "Benchmark is inconsistent.";
* ---------------------------------------------------------------------

         ef = E.L("F");
         pi(r) = PAR_PI.L(r);

**       Store results benchmark
         CTIL(r) = (1/AA(r))*(cprice2(r)/gamma(r))**gamma(r)*(w.L(r)/(1-gamma(r)))**(1-gamma(r));
         CDI_hd = n.l("H")*CTIL("H")/(ttcl("F")*n.l("F")*CTIL("F") + n.l("H")*CTIL("H"));

         BCMRES("Survival cutoff",r,ITR_AA) = z_star.L(r);
         BCMRES("Export cutoff",r, ITR_AA) = zx_star.L(r);
         BCMRES("% Traded, domestic mkt",r,ITR_AA) = 100*(z_star.L(r)/ SUM(rr$oth(r,rr),zx_star.L(rr)))**k(r);
         BCMRES("% Exporters",r,ITR_AA) = 100*(z_star.L(r)/zx_star.L(r))**k(r);
         BCMRES("Mkt share H firms - domestic","H",ITR_AA) = n.l("H")*(1 - theta.L("H"))/(1 - alpha("H"));
         BCMRES("Mkt share H firms - export","H",ITR_AA) = n.l("H")*(1 - theta_x.L("H"))/(1 - alpha("H"));
         BCMRES("Value Exports, bn EUR",r,ITR_AA)       = M.L(r)*n.l(r)*w.L(r)*((z_star.L(r)/zx_star.L(r))**k(r))*(lambda_x(r)/(1-theta_x.L(r)))*coeff_av(r);
         BCMRES("Value Imports, bn EUR","H",ITR_AA)     = M.L("H")*n.l("F")*w.L("F")*((z_star.L("H")/zx_star.L("F"))**k("H"))*(lambda_x("F")/(1-theta_x.L("F")))*coeff_av("H");
         BCMRES("Value Imports, bn EUR","F",ITR_AA)     = M.L("F")*n.l("H")*w.L("H")*((z_star.L("F")/zx_star.L("H"))**k("F"))*(lambda_x("H")/(1-theta_x.L("H")))*coeff_av("F");
         BCMRES("Aggr. (real) Consumption",r,ITR_AA)    =  1/Ptop.L(r)*(W.L(r)*LS(r) + cprice2(r)*E.L(r));
         BCMRES("Aggr. income, bn EUR",r,ITR_AA)        = W.L(r)*LS(r) + cprice2(r)*E.L(r);
         BCMRES("CO2 Emissions, Gt",r,ITR_AA)           = E.L(r);
         BCMRES("Carbon price",r,ITR_AA)                = cprice2(r);
         BCMRES("Carbon costs share",r,ITR_AA)          = cprice2(r)*E.L(r)/(M.L(r)*n.l(r)*Rtil.L(r));
         BCMRES("Markup - Exported",r,ITR_AA)           = 1/theta_x.L(r);
         BCMRES("Markup - Domestic Traded",r,ITR_AA)    = 1/theta.L(r);
         BCMRES("Markup - Domestic Non-Traded",r,ITR_AA)= n.l(r)/(n.l(r) - 1 + alpha(r));
         BCMRES("Share of active varieties",r,ITR_AA)   = M.L(r);
         BCMRES("Number of firms",r,ITR_AA)             = n.l(r);
         BCMRES("Relative wage, Home","H",ITR_AA)       = 1/W.L("F");
         BCMRES("Carbon leakage rate, %","H",ITR_AA)    = na;
         BCMRES("Pass-through rate, Home","H",ITR_AA)   = CDI_hd;
         BCMRES("Pass-through rate, Export","H",ITR_AA) = ttcl("H")*n.l("H")*CTIL("H")/(n.l("F")*CTIL("F") + ttcl("H")*n.l("H")*CTIL("H"));
         BCMRES("Welfare",r,ITR_AA)                     = BCMRES("Aggr. (real) Consumption",r,ITR_AA) - eta_d*SUM(rr, E.L(rr))**2;
         BCMRES("Price level",r,ITR_AA)                 = Ptop.L(r);
         BCMRES("Profits on DG revenues",r,ITR_AA) = (W.L(r)*PAR_FE.L(r)*(1-M.L(r))*n.l(r))/(M.L(r)*n.l(r)*Rtil.L(r));
         BCMRES("Par - cp0",r,ITR_AA) = PAR_cp0.L(r);
         BCMRES("Par - fe",r,ITR_AA) = PAR_FE.L(r);
         BCMRES("Par - lx",r,ITR_AA) = PAR_lx.L(r);
         BCMRES("Par - AA",r,ITR_AA) = AA(r);
         BCMRES("Par - LS",r,ITR_AA) = Ls(r);
         BCMRES("Par - Pi",r,ITR_AA) = PAR_PI.L(r);
         BCMRES("Beta",r,ITR_AA) = beta.L(r);
         BCMRES("Profits",r,ITR_AA) = (W.L(r)*PAR_FE.L(r)*(1-M.L(r))*n.l(r));
         BCMRES("LH_SH",r,ITR_AA) = LH_SH.L(r);
         BCMRES("Emissions av firm",r,ITR_AA)   = E.l(r)/(n.l(r)*M.l(r));
         BCMRES("Emissions av firm f",r,ITR_AA) = m0*gamma(r)*theta.l(r)*w.l(r)*par_l.l(r)*(coeff_av(r)**((1-alpha(r))/alpha(r)))**(sigma(r)-1)/(cprice2(r)*(1-theta.l(r)));
         BCMRES("Emissions av firm x",r,ITR_AA) = m0*gamma(r)*theta_x.l(r)*w.l(r)*par_lx.l(r)*(coeff_av(r)**((1-alpha(r))/alpha(r)))**(sigma(r)-1)/(cprice2(r)*(1-theta_x.l(r)));
         BCMRES("Emissions av firm n",r,ITR_AA) = (BCMRES("Emissions av firm",r,ITR_AA) - omega_out(r)*BCMRES("Emissions av firm f",r,ITR_AA) - {(z_star.l(r)/zx_star.l(r))**k(r)}*BCMRES("Emissions av firm x",r,ITR_AA))/(1-omega_out(r));

         DCMP(ITR_AA,r,"T1") = (theta_a_check(r)/(1-theta_a_check(r)))*PAR_l.L(r)*ratio1_out(r)/z_star.L(r)**(alpha(r)/(1-alpha(r)));
         DCMP(ITR_AA,r,"T2") = (theta.L(r)/(1-theta.L(r)))*PAR_l.L(r)*coeff_av(r);
         DCMP(ITR_AA,r,"T3") = (theta_x.L(r)/(1-theta_x.L(r)))*PAR_lx.L(r)*coeff_av(r);
         DCMP(ITR_AA,r,"1_om") = 1 - omega_out(r);
         DCMP(ITR_AA,r,"om") = omega_out(r);
         DCMP(ITR_AA,r,"chi") = (z_star.L(r)/zx_star.L(r))**k(r);

**       OPTIMAL POLICY
         CTAX.UP = inf;
         CTAX.lo = -inf;
         CTAX.L = scc - PAR_cp0.L("H");

         E.UP(r) = inf;
         E.lo(r) = 0;
         E.L(r) = E0(r);

         beta.UP(r) = 1;
         beta.lo(r) = 0;
         beta.L(r) = shH(r);

**       keep cp0 fixed to benchmark
         PAR_cp0.fx(r) = PAR_cp0.L(r);
         PAR_PI.fx(r) = pi(r);

* ---------------------------------------------------------------------
         SOLVE ASYM_OPT MINIMIZING OBJ_VAR USING NLP;
         abort$(ASYM_OPT.objval > vtol) "Benchmark is inconsistent.";
* ---------------------------------------------------------------------

**       Store results benchmark
         CTIL(r) = (1/AA(r))*(cprice2(r)/gamma(r))**gamma(r)*(w.L(r)/(1-gamma(r)))**(1-gamma(r));
         CDI_hd = n.l("H")*CTIL("H")/(ttcl("F")*n.l("F")*CTIL("F") + n.l("H")*CTIL("H"));

**       Store results simulation of unilateral tax
         OPTRES("Survival cutoff",r,ITR_AA) = z_star.L(r);
         OPTRES("Export cutoff",r, ITR_AA) = zx_star.L(r);
         OPTRES("% Traded, domestic mkt",r,ITR_AA) = 100*(z_star.L(r)/ SUM(rr$oth(r,rr),zx_star.L(rr)))**k(r);
         OPTRES("% Exporters",r,ITR_AA) = 100*(z_star.L(r)/zx_star.L(r))**k(r);
         OPTRES("Mkt share H firms - domestic","H",ITR_AA) = n.l("H")*(1 - theta.L("H"))/(1 - alpha("H"));
         OPTRES("Mkt share H firms - export","H",ITR_AA) = n.l("H")*(1 - theta_x.L("H"))/(1 - alpha("H"));
         OPTRES("Value Exports, bn EUR",r,ITR_AA) = M.L(r)*n.l(r)*w.L(r)*((z_star.L(r)/zx_star.L(r))**k(r))*(lambda_x(r)/(1-theta_x.L(r)))*coeff_av(r);
         OPTRES("Value Imports, bn EUR","H",ITR_AA) = M.L("H")*n.l("F")*w.L("F")*((z_star.L("H")/zx_star.L("F"))**k("H"))*(lambda_x("F")/(1-theta_x.L("F")))*coeff_av("H");
         OPTRES("Value Imports, bn EUR","F",ITR_AA) = M.L("F")*n.l("H")*w.L("H")*((z_star.L("F")/zx_star.L("H"))**k("F"))*(lambda_x("H")/(1-theta_x.L("H")))*coeff_av("F");
         OPTRES("Aggr. (real) Consumption",r,ITR_AA) =  1/Ptop.L(r)*(W.L(r)*LS(r) + cprice2(r)*E.L(r));
         OPTRES("Aggr. income, bn EUR",r,ITR_AA) = W.L(r)*LS(r) + cprice2(r)*E.L(r);
         OPTRES("CO2 Emissions, Gt",r,ITR_AA) = E.L(r);
         OPTRES("Carbon price",r,ITR_AA) = cprice2(r);
         OPTRES("Carbon costs share",r,ITR_AA) = OPTRES("Carbon price",r,ITR_AA)*E.L(r)/(M.L(r)*n.l(r)*Rtil.L(r));
         OPTRES("Markup - Exported",r,ITR_AA) = 1/theta_x.L(r);
         OPTRES("Markup - Domestic Traded",r,ITR_AA) = 1/theta.L(r);
         OPTRES("Markup - Domestic Non-Traded",r,ITR_AA) = n.l(r)/(n.l(r) - 1 + alpha(r));
         OPTRES("Share of active varieties",r,ITR_AA)   = M.L(r);
         OPTRES("Number of firms",r,ITR_AA)             = n.l(r);
         OPTRES("Relative wage, Home","H",ITR_AA) = 1/W.L("F");
         OPTRES("Carbon leakage rate, %","H",ITR_AA) = 100*(E.L("F") - E0("F"))/(E0("H") - E.L("H"));
         OPTRES("Welfare",r,ITR_AA) = OPTRES("Aggr. (real) Consumption",r,ITR_AA) - eta_d*SUM(rr, E.L(rr))**2;
         OPTRES("Price level",r,ITR_AA) = Ptop.L(r);
         OPTRES("Profits on DG revenues",r,ITR_AA) = (W.L(r)*PAR_FE.L(r)*(1-M.L(r))*n.l(r))/(M.L(r)*n.l(r)*Rtil.L(r));
         OPTRES("Par - cp0",r,ITR_AA) = PAR_cp0.L(r);
         OPTRES("Par - fe",r,ITR_AA) = PAR_FE.L(r);
         OPTRES("Par - lx",r,ITR_AA) = PAR_lx.L(r);
         OPTRES("Par - AA",r,ITR_AA) = AA(r);
         OPTRES("Par - Pi",r,ITR_AA) = PAR_PI.L(r);
         OPTRES("Beta",r,ITR_AA) = beta.L(r);
         OPTRES("Profits",r,ITR_AA) = (W.L(r)*PAR_FE.L(r)*(1-M.L(r))*n.l(r));
         OPTRES("LH_SH",r,ITR_AA) = LH_SH.L(r);
         OPTRES("Emissions av firm",r,ITR_AA)   = E.l(r)/(n.l(r)*M.l(r));
         OPTRES("Emissions av firm f",r,ITR_AA) = m0*gamma(r)*theta.l(r)*w.l(r)*par_l.l(r)*(coeff_av(r)**((1-alpha(r))/alpha(r)))**(sigma(r)-1)/(cprice2(r)*(1-theta.l(r)));
         OPTRES("Emissions av firm x",r,ITR_AA) = m0*gamma(r)*theta_x.l(r)*w.l(r)*par_lx.l(r)*(coeff_av(r)**((1-alpha(r))/alpha(r)))**(sigma(r)-1)/(cprice2(r)*(1-theta_x.l(r)));
         OPTRES("Emissions av firm n",r,ITR_AA) = (OPTRES("Emissions av firm",r,ITR_AA) - omega_out(r)*OPTRES("Emissions av firm f",r,ITR_AA) - {(z_star.l(r)/zx_star.l(r))**k(r)}*OPTRES("Emissions av firm x",r,ITR_AA))/(1-omega_out(r));


         DCMP(ITR_AA,r,"T1_OPT") = (theta_a_check(r)/(1-theta_a_check(r)))*PAR_l.L(r)*ratio1_out(r)/z_star.L(r)**(alpha(r)/(1-alpha(r)));
         DCMP(ITR_AA,r,"T2_OPT") = (theta.L(r)/(1-theta.L(r)))*PAR_l.L(r)*coeff_av(r);
         DCMP(ITR_AA,r,"T3_OPT") = (theta_x.L(r)/(1-theta_x.L(r)))*PAR_lx.L(r)*coeff_av(r);
         DCMP(ITR_AA,r,"1_om_OPT") = 1 - omega_out(r);
         DCMP(ITR_AA,r,"om_OPT") = omega_out(r);
         DCMP(ITR_AA,r,"chi_OPT") = (z_star.L(r)/zx_star.L(r))**k(r);


* ===================================================================
* CBAM SCENARIO
* Switch on CBAM to be paied by Foreign
flgcbam("F") = 1;

* ---------------------------------------------------------------------
         SOLVE ASYM_OPT MINIMIZING OBJ_VAR USING NLP;
         abort$(ASYM_OPT.objval > vtol) "Benchmark is inconsistent.";
* ---------------------------------------------------------------------

**       Store results benchmark
         CTIL(r) = (1/AA(r))*(cprice2(r)/gamma(r))**gamma(r)*(w.L(r)/(1-gamma(r)))**(1-gamma(r));
         CDI_hd = n.l("H")*CTIL("H")/(ttcl("F")*n.l("F")*CTIL("F") + n.l("H")*CTIL("H"));

         CBAMRES("Survival cutoff",r,ITR_AA) = z_star.L(r);
         CBAMRES("Export cutoff",r, ITR_AA) = zx_star.L(r);
         CBAMRES("% Traded, domestic mkt",r,ITR_AA) = 100*(z_star.L(r)/ SUM(rr$oth(r,rr),zx_star.L(rr)))**k(r);
         CBAMRES("% Exporters",r,ITR_AA) = 100*(z_star.L(r)/zx_star.L(r))**k(r);
         CBAMRES("Mkt share H firms - domestic","H",ITR_AA) = n.l("H")*(1 - theta.L("H"))/(1 - alpha("H"));
         CBAMRES("Mkt share H firms - export","H",ITR_AA) = n.l("H")*(1 - theta_x.L("H"))/(1 - alpha("H"));
         CBAMRES("Value Exports, bn EUR",r,ITR_AA) = M.L(r)*n.l(r)*w.L(r)*((z_star.L(r)/zx_star.L(r))**k(r))*(lambda_x(r)/(1-theta_x.L(r)))*coeff_av(r);
         CBAMRES("Value Imports, bn EUR","H",ITR_AA) = M.L("H")*n.l("F")*w.L("F")*((z_star.L("H")/zx_star.L("F"))**k("H"))*(lambda_x("F")/(1-theta_x.L("F")))*coeff_av("H");
         CBAMRES("Value Imports, bn EUR","F",ITR_AA) = M.L("F")*n.l("H")*w.L("H")*((z_star.L("F")/zx_star.L("H"))**k("F"))*(lambda_x("H")/(1-theta_x.L("H")))*coeff_av("F");
         CBAMRES("Aggr. (real) Consumption",r,ITR_AA) =  1/Ptop.L(r)*(W.L(r)*LS(r) + cprice2(r)*E.L(r));
         CBAMRES("Aggr. income, bn EUR",r,ITR_AA) = W.L(r)*LS(r) + cprice2(r)*E.L(r);
         CBAMRES("CO2 Emissions, Gt",r,ITR_AA) = E.L(r);
         CBAMRES("Carbon price",r,ITR_AA) = cprice2(r);
         CBAMRES("Carbon costs share",r,ITR_AA) = CBAMRES("Carbon price",r,ITR_AA)*E.L(r)/(M.L(r)*n.l(r)*Rtil.L(r));
         CBAMRES("Markup - Exported",r,ITR_AA) = 1/theta_x.L(r);
         CBAMRES("Markup - Domestic Traded",r,ITR_AA) = 1/theta.L(r);
         CBAMRES("Markup - Domestic Non-Traded",r,ITR_AA) = n.l(r)/(n.l(r) - 1 + alpha(r));
         CBAMRES("Share of active varieties",r,ITR_AA)   = M.L(r);
         CBAMRES("Number of firms",r,ITR_AA)             = n.l(r);
         CBAMRES("Relative wage, Home","H",ITR_AA) = 1/W.L("F");
         CBAMRES("Carbon leakage rate, %","H",ITR_AA) = 100*(E.L("F") - E0("F"))/(E0("H") - E.L("H"));
         CBAMRES("Welfare",r,ITR_AA) = CBAMRES("Aggr. (real) Consumption",r,ITR_AA) - eta_d*SUM(rr, E.L(rr))**2;
         CBAMRES("Price level",r,ITR_AA) = Ptop.L(r);
         CBAMRES("Profits on DG revenues",r,ITR_AA) = (W.L(r)*PAR_FE.L(r)*(1-M.L(r))*n.l(r))/(M.L(r)*n.l(r)*Rtil.L(r));
         CBAMRES("Par - cp0",r,ITR_AA) = PAR_cp0.L(r);
         CBAMRES("Par - fe",r,ITR_AA) = PAR_FE.L(r);
         CBAMRES("Par - lx",r,ITR_AA) = PAR_lx.L(r);
         CBAMRES("Par - AA",r,ITR_AA) = AA(r);
         CBAMRES("Par - Pi",r,ITR_AA) = PAR_PI.L(r);
         CBAMRES("Beta",r,ITR_AA) = beta.L(r);
         CBAMRES("Profits",r,ITR_AA) = (W.L(r)*PAR_FE.L(r)*(1-M.L(r))*n.l(r));
         CBAMRES("LH_SH",r,ITR_AA) = LH_SH.L(r);
         CBAMRES("Emissions av firm",r,ITR_AA)   = E.l(r)/(n.l(r)*M.l(r));
         CBAMRES("Emissions av firm f",r,ITR_AA) = m0*gamma(r)*theta.l(r)*w.l(r)*par_l.l(r)*(coeff_av(r)**((1-alpha(r))/alpha(r)))**(sigma(r)-1)/(cprice2(r)*(1-theta.l(r)));
         CBAMRES("Emissions av firm x",r,ITR_AA) = m0*gamma(r)*theta_x.l(r)*w.l(r)*par_lx.l(r)*(coeff_av(r)**((1-alpha(r))/alpha(r)))**(sigma(r)-1)/(cprice2(r)*(1-theta_x.l(r)));
         CBAMRES("Emissions av firm n",r,ITR_AA) = (CBAMRES("Emissions av firm",r,ITR_AA) - omega_out(r)*CBAMRES("Emissions av firm f",r,ITR_AA) - {(z_star.l(r)/zx_star.l(r))**k(r)}*CBAMRES("Emissions av firm x",r,ITR_AA))/(1-omega_out(r));

         DCMP(ITR_AA,r,"T1_CBAM")   = (theta_a_check(r)/(1-theta_a_check(r)))*PAR_l.L(r)*ratio1_out(r)/z_star.L(r)**(alpha(r)/(1-alpha(r)));
         DCMP(ITR_AA,r,"T2_CBAM")   = (theta.L(r)/(1-theta.L(r)))*PAR_l.L(r)*coeff_av(r);
         DCMP(ITR_AA,r,"T3_CBAM")   = (theta_x.L(r)/(1-theta_x.L(r)))*PAR_lx.L(r)*coeff_av(r);
         DCMP(ITR_AA,r,"1_om_CBAM") = 1 - omega_out(r);
         DCMP(ITR_AA,r,"om_CBAM")   = omega_out(r);
         DCMP(ITR_AA,r,"chi_CBAM")  = (z_star.L(r)/zx_star.L(r))**k(r);


* Turn CBAM off again.
flgcbam(r) = 0;

);


CHANGE_OPT("CO2 Emissions, %","Pigou",r,ITR_AA) = 100*(OPTRES("CO2 Emissions, Gt",r,ITR_AA) - BCMRES("CO2 Emissions, Gt",r,ITR_AA))/BCMRES("CO2 Emissions, Gt",r,ITR_AA);
CHANGE_OPT("Share of active varieties, %","Pigou",r,ITR_AA) = 100*(OPTRES("Share of active varieties",r,ITR_AA) - BCMRES("Share of active varieties",r,ITR_AA))/BCMRES("Share of active varieties",r,ITR_AA);
CHANGE_OPT("Number of firms, %","Pigou",r,ITR_AA) = 100*(OPTRES("Number of firms",r,ITR_AA) - BCMRES("Number of firms",r,ITR_AA))/BCMRES("Number of firms",r,ITR_AA);
CHANGE_OPT("Emissions av firm n, %","Pigou",r,ITR_AA) = 100*(OPTRES("Emissions av firm n",r,ITR_AA) - BCMRES("Emissions av firm n",r,ITR_AA))/BCMRES("Emissions av firm n",r,ITR_AA);
CHANGE_OPT("Emissions av firm f, %","Pigou",r,ITR_AA) = 100*(OPTRES("Emissions av firm f",r,ITR_AA) - BCMRES("Emissions av firm f",r,ITR_AA))/BCMRES("Emissions av firm f",r,ITR_AA);
CHANGE_OPT("Emissions av firm x, %","Pigou",r,ITR_AA) = 100*(OPTRES("Emissions av firm x",r,ITR_AA) - BCMRES("Emissions av firm x",r,ITR_AA))/BCMRES("Emissions av firm x",r,ITR_AA);
CHANGE_OPT("omega, %","Pigou",r,ITR_AA) = 100*(DCMP(ITR_AA,r,"om_OPT") - DCMP(ITR_AA,r,"om"))/DCMP(ITR_AA,r,"om");
CHANGE_OPT("chi, %","Pigou",r,ITR_AA) = 100*(DCMP(ITR_AA,r,"chi_OPT") - DCMP(ITR_AA,r,"chi"))/DCMP(ITR_AA,r,"chi");
CHANGE_OPT("Markup - n, %","Pigou",r,ITR_AA) = 100*(OPTRES("Markup - Domestic Non-Traded",r,ITR_AA) - BCMRES("Markup - Domestic Non-Traded",r,ITR_AA))/BCMRES("Markup - Domestic Non-Traded",r,ITR_AA);
CHANGE_OPT("Markup - f, %","Pigou",r,ITR_AA) = 100*(OPTRES("Markup - Domestic Traded",r,ITR_AA) - BCMRES("Markup - Domestic Traded",r,ITR_AA))/BCMRES("Markup - Domestic Traded",r,ITR_AA);
CHANGE_OPT("Markup - x, %","Pigou",r,ITR_AA) = 100*(OPTRES("Markup - Exported",r,ITR_AA) - BCMRES("Markup - Exported",r,ITR_AA))/BCMRES("Markup - Exported",r,ITR_AA);
CHANGE_OPT("Value Exports, %","Pigou",r,ITR_AA) = 100*(OPTRES("Value Exports, bn EUR",r,ITR_AA) - BCMRES("Value Exports, bn EUR",r,ITR_AA))/BCMRES("Value Exports, bn EUR",r,ITR_AA);
CHANGE_OPT("Welfare, %","Pigou",r,ITR_AA) = 100*(OPTRES("Welfare",r,ITR_AA) - BCMRES("Welfare",r,ITR_AA))/BCMRES("Welfare",r,ITR_AA);


CHANGE_CBAM("CO2 Emissions, %","Pigou",r,ITR_AA) = 100*(CBAMRES("CO2 Emissions, Gt",r,ITR_AA) - BCMRES("CO2 Emissions, Gt",r,ITR_AA))/BCMRES("CO2 Emissions, Gt",r,ITR_AA);
CHANGE_CBAM("Share of active varieties, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Share of active varieties",r,ITR_AA) - BCMRES("Share of active varieties",r,ITR_AA))/BCMRES("Share of active varieties",r,ITR_AA);
CHANGE_CBAM("Number of firms, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Number of firms",r,ITR_AA) - BCMRES("Number of firms",r,ITR_AA))/BCMRES("Number of firms",r,ITR_AA);
CHANGE_CBAM("Emissions av firm n, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Emissions av firm n",r,ITR_AA) - BCMRES("Emissions av firm n",r,ITR_AA))/BCMRES("Emissions av firm n",r,ITR_AA);
CHANGE_CBAM("Emissions av firm f, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Emissions av firm f",r,ITR_AA) - BCMRES("Emissions av firm f",r,ITR_AA))/BCMRES("Emissions av firm f",r,ITR_AA);
CHANGE_CBAM("Emissions av firm x, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Emissions av firm x",r,ITR_AA) - BCMRES("Emissions av firm x",r,ITR_AA))/BCMRES("Emissions av firm x",r,ITR_AA);
CHANGE_CBAM("omega, %","Pigou",r,ITR_AA) = 100*(DCMP(ITR_AA,r,"om_CBAM") - DCMP(ITR_AA,r,"om"))/DCMP(ITR_AA,r,"om");
CHANGE_CBAM("chi, %","Pigou",r,ITR_AA) = 100*(DCMP(ITR_AA,r,"chi_CBAM") - DCMP(ITR_AA,r,"chi"))/DCMP(ITR_AA,r,"chi");
CHANGE_CBAM("Markup - n, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Markup - Domestic Non-Traded",r,ITR_AA) - BCMRES("Markup - Domestic Non-Traded",r,ITR_AA))/BCMRES("Markup - Domestic Non-Traded",r,ITR_AA);
CHANGE_CBAM("Markup - f, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Markup - Domestic Traded",r,ITR_AA) - BCMRES("Markup - Domestic Traded",r,ITR_AA))/BCMRES("Markup - Domestic Traded",r,ITR_AA);
CHANGE_CBAM("Markup - x, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Markup - Exported",r,ITR_AA) - BCMRES("Markup - Exported",r,ITR_AA))/BCMRES("Markup - Exported",r,ITR_AA);
CHANGE_CBAM("Value Exports, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Value Exports, bn EUR",r,ITR_AA) - BCMRES("Value Exports, bn EUR",r,ITR_AA))/BCMRES("Value Exports, bn EUR",r,ITR_AA);
CHANGE_CBAM("Welfare, %","Pigou",r,ITR_AA) = 100*(CBAMRES("Welfare",r,ITR_AA) - BCMRES("Welfare",r,ITR_AA))/BCMRES("Welfare",r,ITR_AA);


DECOMP("dt",r)             = OPTRES("Carbon price",r,"ITR_AA4")/BCMRES("Carbon price",r,"ITR_AA4")  - 1;
DECOMP("weight n firms",r) = DCMP("ITR_AA4",r,"1_om") * BCMRES("Emissions av firm n",r,"ITR_AA4") / BCMRES("Emissions av firm",r,"ITR_AA4");
DECOMP("dtheta_n",r)       = (BCMRES("Markup - Domestic Non-Traded",r,"ITR_AA4")/OPTRES("Markup - Domestic Non-Traded",r,"ITR_AA4") -1 )  / (1-1/BCMRES("Markup - Domestic Non-Traded",r,"ITR_AA4"));
DECOMP("weight f firms",r) = DCMP("ITR_AA4",r,"om") * BCMRES("Emissions av firm f",r,"ITR_AA4") / BCMRES("Emissions av firm",r,"ITR_AA4");
DECOMP("dtheta_f",r)       = (BCMRES("Markup - Domestic Traded",r,"ITR_AA4")/OPTRES("Markup - Domestic Traded",r,"ITR_AA4") -1 ) / (1-1/BCMRES("Markup - Domestic Traded",r,"ITR_AA4"));
DECOMP("weight x firms",r) = DCMP("ITR_AA4",r,"chi") * BCMRES("Emissions av firm x",r,"ITR_AA4") / BCMRES("Emissions av firm",r,"ITR_AA4");
DECOMP("dtheta_x",r)       = (BCMRES("Markup - Exported",r,"ITR_AA4")/OPTRES("Markup - Exported",r,"ITR_AA4") -1 ) / (1-1/BCMRES("Markup - Exported",r,"ITR_AA4"));
DECOMP("domega_n",r)       = (CHANGE_OPT("omega, %","Pigou",r,"ITR_AA4")/100) * DCMP("ITR_AA4",r,"om")/DCMP("ITR_AA4",r,"1_om");
DECOMP("domega_f",r)       = (CHANGE_OPT("omega, %","Pigou",r,"ITR_AA4")/100);
DECOMP("dchi_x",r)         = (CHANGE_OPT("chi, %","Pigou",r,"ITR_AA4")/100);


display CHANGE_OPT, DECOMP

* ---------------------------------------------------------------------
*** Figures and tables in the paper
* Figure 2 shows results stored in DECOMP
* Tables 3 and 4 show results in CHANGE_OPT, all iterations
* Tables 5 and 6 show results in CHANGE_CBAM, iteration ITR_AA4

* ---------------------------------------------------------------------
DISPLAY BCMRES, OPTRES, DCMP;

IF (PRINT = 1,


execute_unload "CHANGE_NOV21.gdx", CHANGE_OPT, CHANGE_CBAM
execute 'gdxxrw.exe CHANGE_NOV21.gdx par=CHANGE_OPT cdim=3 rng=CHN_OPT!B2';
execute 'gdxxrw.exe CHANGE_NOV21.gdx par=CHANGE_CBAM cdim=3 rng=CHN_CBAM!B2';

execute_unload "ALLRES_NOV21.gdx", BCMRES OPTRES CBAMRES
execute 'gdxxrw.exe ALLRES_NOV21.gdx par=BCMRES cdim=2 rng=BCM!A1';
execute 'gdxxrw.exe ALLRES_NOV21.gdx par=OPTRES cdim=2 rng=OPT!A1';
execute 'gdxxrw.exe ALLRES_NOV21.gdx par=CBAMRES cdim=2 rng=CBAM!A1';


execute_unload "DCMP_NOV21.gdx", DECOMP
execute 'gdxxrw.exe DCMP_NOV21.gdx par=DECOMP cdim=1 rng=data!A1';

execute_unload "DECOMP_NOV21.gdx", DECOMP
execute 'gdxxrw.exe DECOMP_NOV21.gdx par=DECOMP cdim=1 rng=data!A1';


);
* ---------------------------------------------------------------------







*$OFFTEXT
